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NMR SIMULATOR 


DECUS Program Library Write-up DECUS No. 8-194 

ABSTRACT 

NMR Simulation is designed to calculate the theoretical spectrum of compounds containing 
hydrogen, fluorine, carbon-13 and other nuclei of spin 1/2. The calculated theoretical spectrum 
is displayed on an oscilloscope. 

Options for punched and typewritten output, change in x-axis offset (sweep offset) and spectrum 
resolution are available. Chemical shift and coupling constant parameters may be varied succes¬ 
sively until the display spectrum matches that obtained experimentally. Redisplay of a "library" 
of theoretical spectra is possible by retaining punched output tapes. 

REQUIREMENTS 

Equipment 


A PD P—8/1 equipped with high-speed reader and punch. Teletype, 8K of memory, and AX08 
are required. (Without AX08 and/or high-speed reader and punch see Appendix A.) 

Memory 

Field J0-program uses 0 - 4577 and 5400 - 7577. Field 1-data is stored in location 20 - 7245. 
Loading Procedure 


Load the program into field 0 with the standard binary loader. 
Operating Procedure 


Set SR - 000000 press STOP, press LOAD ADDRESS, and press START. 

The computer then begins to ask for information. The Teletype questions along with appropriate 
operator action are shown in Table I. (See Appendix B for typing error recovery.) 
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TABLE I 

ORDER OF ENTERING DATA 


Section Teletype Output 

A "COMMENTS" 


B "THE NUMBER OF SPINS IS" 


C "PUNCHED OUTPUT?" 


D "SWEEP OFFSET:" 


E "SWEEP WIDTH:" 


F "CHEMICAL SHIFTS (CPS)" 


G "COUPLING CONSTANTS (CPS)" 


Operator Action 


Type any amount of commentary such as 
the problem name etc. Exit to the rest 
of the program is achieved by pressing 
ALT MODE on the Teletype. 


Enter the number of spins. Values from 
2-6 will be accepted. With other 
values, "? " will be typed and the question 
asked again. 


Type Y if punched output is desired, or 
N if it is not. If Y is typed, turn on the 
high-speed punch. 


Enter the desired sweep offset (in cps). 
TYPE RETURN. This is the extent in cps 
to which the display is shifted to the right. 

Type the number of cps to be displayed on 
the scope followed by a RETURN. Smaller 
width gives greater resolution. (Entering a 
0 leads to division by 0 and one must return 
to SWEEP OFFSET (See below).) 


Enter the chemical shifts (in cps) in the 
order v (1), v (2),...., v (n) 
following each with a RETURN. (See 
appendix C.) 


Enter the coupling constants (in cps) in the 


° rder J l,2' J l,3' 




J 0 ,J 0 „...,J , following each 

J,4 n-l,n 

with a RETURN. (See appendix C.) 


The portion of the spectrum specified by the sweep width and sweep offset will now be displayed on 
the scope after the calculation is complete. (See appendix C.) It will be noted that, at regular 
intervals along the baseline, bright spots appear. These are calibration marks. The number of cps 
between two calibration marks is 1/25 of the sweep width. Thus, at a sweep width of 250 cps, there 
are 10 cps between calibration marks, but, at 50 cps sweep width, there are only 2 cps between 
calibration marks. 
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While the spectrum is being displayed, the operator can type any amount of commentary. When 
editing of the displayed spectrum is desired, typing one of the mnemonics in the table below followed 
by pressing ALT MODE on the Teletype will achieve exit to the portion of the program described in 
Table II. 


TABLE II 


COMMANDS AVAILABLE 


OPTION 

# 


1 



2 

3 


4 


5 


6 


MNEMONIC (to be typed 
by operator) 

NEW DATA 

START OVER 

NEW COUPLING CONSTANTS 
OFFSET 

DO AGAIN 

LISTING 


EXIT TO (The letters refered to program 
sections in Table I) 


Section F. Enter new chemical shifts and 
coupling constants. 

Section A. Begin a new problem. 

Section G. Enter new coupling constants. 

See below. Allows display of a different 
portion of the spectrum. 

Section D. Enter new sweep offset, 
sweep width, chemical shifts, and coupling 
constants. 

See below. Allows a listing of transition 
energies and intensities on Teletype. (A 
punched tape is required for this option.) 


• Option 4: If this option is selected only parts D and E of Table I are repeated. After entering a 
new sweep offset and sweep width, the Teletype will print "TAPE?". If a punched tape has been 
prepared, type Y. If not, type N. If N is selected, the spectrum will be recalculated and displayed. 
If Y has been typed, the Teletype will write "LOAD AND CONT.". Place the tape in the high-speed 
reader and press CONTINUE. After the tape has been read, the offset spectrum will be displayed. 

If it is desired to display a previously punched tape, the program may be started at 000200. This 
leads directly into the offset routine. When starting at this address, only options 2, 4, or 6 should 
be selected when the display spectrum appears. 


Option 6: If this option is selected, "MINIMUM INTENSITY:" will be typed. Enter the minimum 
intensity to be listed. (Transitions with intensities lower than this will not be listed.) Type RETURN. 
The Teletype will then print "LOAD AND CONT.". Place the punched tape in the high-speed 
reader and press CONTINUE. "ENERGY INTENSITY" will then be typed, and the transition energies 
and intensities will be listed. When the listing is done, the spectrum will be redisplayed on the scope 
at half its original intensity. If a listing of a previously punched tape is desired, the program may 
be started at 003400. This leads directly to the listing routine. Again, when starting at this address 
only options 2, 4, or 6 should be selected when the spectrum appears. When the program is started 
at this address, there is no way of predicting what will be displayed on the scope. 


3 






Operating Suggestions 
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The punched output option (Section C Table I) is designed to provide a permanent record of a 
calculated spectrum, for display at any future time, and to provide listings of transition energies 
and intensities. Some time is saved in the five and six spin systems by having a tape punched and 
then using this for inspection with several offsets and different resolutions (i.e. Sweep Widths). 
However, the original spectrum calculation is rapid enough so that this is not saving for smaller 
systems. Thus it is probably best to select the N option here unless a large system (at least 5 spins) 
is being calculated or unless a permanent record is desired. If Y is selected, tapes will be punched 
until the N option is selected either by restarting the program or by using option 2 of Table II. 

A useful starting point for many 60 megacycle nmr (hydrogen) spectra of organic molecules is a 
sweep offset of 0 cps and a sweep width of 500 cps. For best results, however, all peaks should 
eventually be examined at a considerably smaller sweep width, since rounding errors may give 
misleading scope displays with large sweep widths. It should be noted that increasing the sweep offset 
moves peaks to the right, and decreasing the sweep offset moves peaks to the left. Once the peaks 
are observed, the proper sweep width and sweep offset can be determined from the calibration marks. 

As a starting point, 0.02 is a good intensity to select as the minimum intensity for option 6 Table II. 
Typing 0 will allow listing of all transitions, and for small systems this is wise. This is also the 
approximate minimum intensity that will be displayed on the scope. 

Appendix D contains some typical input output from this program. 

NMR Background 





The discussion below assumes that the nuclei being studied are hydrogen, fluorine, carbon-13 or 
other nuclei of spin 1/2. 

Each proton (or fluorine nucleus) in a molecule can be assigned a spin of either m or P . For a 
system containing n nuclei, there are 2 n spin states possible. Each of these states is called a basic 
product function. The energy of a basic product function can be calculated from equation (1); where 
s; = -1/2 for |3 spin and + 1/2 for a spin; Vj is the frequency of absorption in cps; T;j = + 1/4 if 
i and j have the same spin and -1/4 if i and j have opposite spin; and Jj- is the coupling constant 
between i and j in cps. 

(1) H uu= E (siV |+ E TjjJjj) 

' = 1 j>i 

The energy of interaction between two basic product functions is 0 unless both have the same number 
of a and p spins. When this condition is fulfilled, the interaction energy is calculated from equation 
(2) where U = 1 if the basic product functions differ only in the interchange of spins i and j. (J;j as 
defined above.) Otherwise, U =0. 

(2) H uv = 1/2UJ,, 

Thus, to calculate an mnr spectrum, one contructs all of the possible basic product functions and sorts 
them into groups, each member of a group containing the same number of a and (3 spins. The energy 
of each of them is determined and used as the diagnal element in a square matrix. The off-diagonal 
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elements of this matrix (H matrix) are energies of interaction between the members of the group. 
This matrix is diagonalized and the trace of the matrix then contains the energies of the final spin 
functions. These final spin eigenfunctions are the columns of the matrix (U matrix) required to 
diagonalize the H matrix. Each element in a spin function represents the contribution of a basic 
product function to that spin function. 


Before continuing, it is necessary to define the F z value of a spin function as shown in equation 
(3) where the s; is as previously defined. 

n 


(3) 


F = 
z 


i= 1 


Now the peaks on a recorder nmr spectrum are transitions from one spin function to another. These 
transition energies and their intensities can be calculated as outlined below. Transitions are allowed 
between spin functions whose F z values differ only by 1. If this holds, the energy of the transition 
is the difference in energy of the final spin eigenfunctions, and the intensity of the transition is given 
by equation (4) where is the uth element in one spin function, is the vth element in the other 
spin function, and A = 1 if the basic product functions represented by C (j and differ by one spin. 
Otherwise; A = 0. 


(4) I = ( Z 1 C C ' A) 2 

x/ UV U V ' 

Program Description 


The complete program listings are contained in Appendix E with extensive commentary. The following 
is a brief outline of the program; capital letters below refer to variables. 


I. 


IV. 


V. 

VI. 

VII. 

VIII. 


IX. 


X. 

XI. 


XII. 

XIII. 


Begins with text mode for input of commentary. 

Enter the number of spins and store in N. 

Set "switch" for punched output. 

Set up basic product functions (BPF's). 

A. Calculate NARRY array which contains relative addresses of the BPF's of different 
F value. 

B. Cal culate NUSE array. This array actually contains the BPF's. They are stored in 
groups. Each group having an F z value one greater than the preceeding group. 

I 

Get sweep offset and sweep width. 

Get chemical shifts. 

Get coupling constants. 

Initialize SPEC array. 

A. This array eventually contains the displayed spectrum. 

B. It is intialized to 0 except for calibration points where it is set to 4000g. 

Calculate first H and U matrices. 

A. First H matrix is always a 1 X 1 which does not require diagonalization. 

B. First U matrix is always a 1 X 1 with its only element equal to 1. 

Set FLAG which is the negative of the number of times the loop from XI to XVII must be 
executed. 

Set N2, calculate EN, and UOLD arrays. 

A. Copy the trace of the H matrix into the EN array. 

B. Copy the U matrix into the UOLD array. 

C. Set N2 equal to the size of the H matrix just copied. 

Calculate the size of the next H matrix and put it into N1. 

Determine which transitions will be allowed between the BPF's represented by UOLD and those 
that will be represented by U. This is stored in the array called TABLE. 
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XIV. 


XV. 

XVI. 


XVII. 

XVIII. 


Initialize H and U matrices. 

A. Diagonal H elements are calculated by equation (1). 

B. Off-diagonal H elements are calculated by equation (2). 

C. The U matrix is initialized to all zeros except for ones on the diagonal. 
Diagonalize the H matrix and calculate the U matrix by the Jacobi method. 
Calculate transitions. 

A. Calculate transition energies by forming the difference between the elements 
of the trace of the H matrix and the elements in the EN array. 

B. Calculate transition intensities from equation (4). 

C. Store appropriate values in SPEC array. 

Test FLAG to see if the loop is done. 

A. If loop is not done, go to XI. 

If loop is done, display spectrum. 

If "ALT MODE" is typed, interpret mnemonic and execute option. 


The listings include everything except Floating Point Package II (DEC-08-YQ2A-PB). 

Included, however, is an overlay of Floating Point Package II. These modifications force all 
indirect floating commands to operate on field 1. The field width for the output controller is now 

stored at location 3 in this version. All direct commands of the floating point package still operate 
on field 0. 
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APPENDIX A 

If the AX08 is not available, the table below shows the location of the scope control commands 
and their function. These then may be replaced with their counterpart instructions. 


Location 

Contents 

Function 

4511 

6302 

clear and load X buffer 

4516 

6317 

clear and load Y buffer, and display 

4527 

6317 

clear and load Y buffer, and display 

4541 

6317 

clear and load Y buffer, and display 

4542 

6317 

clear and load Y buffer, and display 

4543 

6317 

clear and load Y buffer, and display 

ligh-speed reader and punch are not available 

the following changes should be made: 

Location 

Contents 

Change To 

0202 

6016 

7000 

024] 

6014 

6032 

0244 

6011 

6031 

0246 

6016 

6036 

3135 

6026 

6046 

3325 

4336 

4134 

3334 

4336 

4134 

3402 

6016 

7000 

3447 

6014 

6032 

3473 

6011 

6031 

3475 

6016 

6036 

4017 

6014 

7000 

4020 

6026 

7000 

4272 

6021 

6041 

4274 

6026 

6046 

4323 

6011 

6031 

4325 

6016 

6036 

4337 

6011 

6031 

4341 

6016 

6034 

4353 

6026 

6046 

4351 

7200 

7602 

4322 

7300 

6032 

4326 

7006 

7106 


When program is started punch should be off. Enter data as described in Table I. After the last 
coupling constant has been entered the computer will halt. Turn punch on, turn Teletype to local, 
press "HERE IS" on Teletype, turn Teletype to line, and press CONTINUE. Simply press CONTINUE 
if punched output is not desired. 
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APPENDIX B 

Location of Typing 

Error 

Recovery 

Section B Table 1 

Any character typed will be interpreted. If an error is 
made, the Teletype will echo the character as it was 
interpreted (last ASCII code digit). Thus, if "D" is typed, 

"4" will be echoed and stored as the number of spins. 

Recovery can be made by restarting program. 

Section C Table 1 

Restart program. 

Section D, E, F, and 

G Table 1, "MINIMUM 
INTENSITY:" of option 

6 Table II 

These use the standard floating point input. If an error is 
made, type rubout, and retype the entire number. An ^ 

illegal character will terminate the input and recovery can 
be made by restarting the program or continuing until the 
display appears and selecting a suitable option from Table II. 

"TAPE?" in option 4 

Table II 

Restart program at 000200, 

Mnemonic during display 

Type RETURN, and retype mnemonic. 
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EXAMPLE I 

COMMENTS 

SAMPLE AB PATTERN, CHEMICAL SHIFTS 265 AND 280CPS, COUPLING 
CONSTANT 5 CPS 

THE NUMBER OF SPINS IS 2 

PUNCHED OUTPUT? Y 

SWEEP OFFSET: 0 

SWEEP WIDTH: 500 

CHEMICAL SHIFTS (CPS) 

265 

280 

COUPLING CONSTANTS (CPS) 

5 

THE DISPLAYED SPECTRUM IS SHOWN HERE: 


FT" I i r 'l i i i i i —i—i—i— rm i i i—i—i— 

NOW USE THE OFFSET OPTION TO ENLARGE THE DISPLAY: 
OFFSET 

SWEEP OFFSET: 250 

SWEEP WIDTH: 50 

TAPE? Y LOAD AND CONT. 

THE DISPLAYED SPECTRUM IS SHOWN BELOW: 


r i i i i i i i 'i i I i — r i i I ii — r T » i i » 
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CN CO ^ lO o 


I 


^ of Spins 
Selected 


APPENDIX C 


# of Chemical ^ of Coupling Execution Time (SEC) 

Shifts Required Constants Required if Tape is Not Punched 


2 

3 

4 

5 

6 


1 

3 

6 

10 

15 


<1 

2-3 

8-10 

30-90 

300-600 
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THE LISTING OPTION WILL TYPE OUT THE ENERGIES AND INTENSITIES OF THE SPECTRUM 
WHICH WAS DISPLAYED: 

LISTING 


MINIMUM INTENSITY: 0.02 
LOAD AND CONT. 


ENERGY 


INTENSITY 

+ 277.905 

+ 

1.316 

+ 262.094 

+ 

.683 

+ 267.094 

+ 

1.316 

+ 282.905 

+ 

.683 


THE SCOPE DISPLAY IS NOW EXACTLY THE SAME AS BEFORE EXCEPT THAT THE INTENSITY HAS 
BEEN REDUCED BY HALF. 




^ 00000000^00 00 ^^ 


EXAMPLE 2 


A COMPLETELY DIFFERENT SYSTEM WILL NOW BE CALCULATED: 
START OVER 


COMMENTS 

ETHYL GROUP, CHEMICAL SHIFTS 90 AND 260CPS, COUPLING CONSTANTS 8 CPS 

THE NUMBER OF SPINS IS 5 

PUNCHED OUTPUT? N 

SWEEP OFFSET: 0 

SWEEP WIDTH: 500 

CHEMICAL SHIFTS (CPS) 

90 

90 

90 

260 

260 

COUPLING CONSTANTS (CPS) 


THE DISPLAYED SPECTRUM IS SHOWN BELOW: 
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Program Listing 


The listing is organized according to the table of contents shown below. 

Starting address 


Contents 

(Octal) 

listing 

Common constants, read and write routine, 
base addresses of arrays 

0 

1 

Routine to offset spectrum 

200 

2 

Calculate NARRAY and NUSE arrays 

400 

4 

Input chemical shifts 

600 

6 

Input coupling constants 

610 

6 

Calculate diagonal Helements 

645 

6 

Calculate off-diagonal Helements 

1000 

9 

This routine is not used 

1130 

10 

Initialize SPEC array 

1151 

10 

Search H array to find elements requiring reduction 

1200 

12 

Calculate position in H array from subscripts 

1273 

13 

Calculate position in U array from subscripts 

1341 

13 

Calculate parameters to convert one H element to 0 

1400 

15 

Alter H and U elements 

1600 

17 

Calculate TABLE array which contains which transitions 
are allowed 

2000 

19 

Get initial values for H and U matrices 

2200 

21 

Copy trace elements of H matrix into EN array and 

U into UOLD array 

2325 

22 

Arrange calculation of transition energy and intensity 

2400 

24 

Organize punching of floating accumulator 

2554 

26 

Routine to normalize and combine adjacent peaks for 
display 

2600 

27 
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Contents 

Starting address 
(Octal) 

Page of 
listing 

Routine to accept commentary and mnemonics 
during display 

3000 

29 

Routine for alphanumeric output 

3057 

29 

Routine to calculate transition energy 

3200 

31 

Routine calculates transition probability and increments 3247 

SPEC array 

31 

Routine to punch contents of accumulator 

3315 

32 

Routine to fix a floating point number 

3345 

32 

Routine to list energies and intensities 

3400 

34 

Main driver program outlined on pp. 5-6 

4000 

37 

Routine to read a number from tape into floating 
accumulator 

4314 

40 

Routine to get sweep offset and sweep width 

4400 

42 

Routine to display spectrum 

4500 

43 

Floating point modification 

— 

44 
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